<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Visualize Elementary Flux Modes from Human model </title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S2 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S3 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S4 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: pre-wrap;  }
.S5 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S6 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S7 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S9 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S12 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: pre;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: pre;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S13 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>Visualize Elementary Flux Modes from Human model </span></h1><h2  class = 'S1'><span>Author(s): Chaitra Sarathy (Maastricht University, The Netherlands)</span></h2><div  class = 'S2'><span style=' font-style: italic;'></span></div><h2  class = 'S1'><span>Reviewer(s): </span></h2><div  class = 'S2'><span style=' font-style: italic;'>Include those who have reviewed the tutorial, and are therefore partly responsible for picking up errors or omissions.</span></div><h2  class = 'S1'><span>INTRODUCTION</span></h2><div  class = 'S2'><span>Elementary Flux Modes (EFMs) are a tool for constraint-based modelling and metabolic network analysis. However, systematic and automated visualization of EFMs, capable of integrating various data types is still a challenge. This tutorial is part of the study by Sarathy et al. (2019, in preparation), "EFMviz: A COBRA Toolbox extension to visualize Elementary Flux Modes in genome-scale metabolic models". </span></div><div  class = 'S2'><span>This workflow can be used to visualize EFMs calculated using any currently available methods. Inputs for the workflow are: (1) the text file containing EFMs and (2) the genome-scale metabolic model used for generating EFMs. Additionally, omics data can be used to select one (or more) EFM(s) for visualization among the vast number of EFMs that usually results from EFM calculation. Sarathy et al. (2019, in preparation) demonstrate two examples of such a selection procedure (a) Yield Analysis and (b) EFM Enrichment with EFMs calculated from genome-scale metabolic models, (1) iAF1260 (E. coli) [1] and (2) Recon 2.2 (Human) [2]. </span></div><div  class = 'S2'><span>Each use case is explained with a tutorial. This tutorial corresponds to the second use case and it describes how to analyze and visualize EFMs from the human</span><span style=' font-style: italic;'> </span><span>model Recon 2.2</span><span style=' font-style: italic;'>. </span><span>The aim is to apply the workflow on a large-scale model and it can be used as a means to study cancer metabolism. EFMs were generated for glucose uptake reaction. For details on how EFMs were generated, please refer the citation section at </span><a href = "https://github.com/macsbio/efmviz"><span>https://github.com/macsbio/efmviz</span></a></div><div  class = 'S2'><span>The steps in the workflow include:</span></div><div  class = 'S2'><span>        </span><span style=' font-weight: bold;'>A. Specifying Inputs</span></div><div  class = 'S2'><span>        i. Import input data: file containing EFMs (.txt), model file (.mat), file containing gene level statistics typically resulting from a differential gene expression analysis (.txt)</span></div><div  class = 'S2'><span>       </span><span style=' font-weight: bold;'> B. EFM Selection</span></div><div  class = 'S2'><span>        i. Filter EFMs containing a reaction of interest</span></div><div  class = 'S2'><span>        ii. Perform EFM Enrichment Analysis</span></div><div  class = 'S2'><span>            a. Generate input files for EFM Enrichment (preprocessing)</span></div><div  class = 'S2'><span>            b. Call R script for Enrichment (EFM_enrichment.R)</span></div><div  class = 'S2'><span>        iii. EFM Selection</span></div><div  class = 'S2'><span>            a. Get indices of enriched EFMs </span></div><div  class = 'S2'><span>            b. Find Subsystem Occurrence in enriched EFMs for selecting an EFM for visualization, for example, EFM with reactions in TCA Cycle</span></div><div  class = 'S2'><span>        </span><span style=' font-weight: bold;'>C. Submodel Creation</span></div><div  class = 'S2'><span>        i. Write submodel (with ubiquitous metabolites)</span></div><div  class = 'S2'><span>        ii. Write submodel (withOUT ubiquitous metabolites)</span></div><div  class = 'S2'><span>        iii. Extract reactions occurring in a certain percentage of EFMs</span></div><div  class = 'S2'><span>        </span><span style=' font-weight: bold;'>D. Visualization</span></div><div  class = 'S2'><span>        Open the SBML file in Cytoscape and map the data used during step B.</span></div><div  class = 'S2'><span>The following sections will walk you through the code for performing the steps A-D, thereby reproducing the main results from Sarathy et al (2019, in preparation) (Use Case 2). For more details on the workflow please refer </span><a href = "https://github.com/macsbio/efmviz"><span>https://github.com/macsbio/efmviz</span></a></div><div  class = 'S2'><span style=' font-weight: bold;'>Input File Formats</span></div><ul  class = 'S3'><li  class = 'S4'><span>model (.mat): same as the model used for EFM generation </span></li><li  class = 'S4'><span>file containing EFMs (.txt): each row corresponds to an EFM and column has indices of reactions active in the EFM. Each row can have varied number of columns. The first column should contain an EFM ID.</span></li><li  class = 'S4'><span>file containing gene level statistics (.txt): Tab separated file containing an expression level (log2 fold change) and p-value for every gene. Columns: &lt;GeneID&gt;&lt;GeneName&gt;&lt;log2FC&gt;&lt;PValue&gt;. Make sure to have the gene IDs in the same format as in the model.</span></li></ul><h2  class = 'S5'><span>PROCEDURE</span></h2><h2  class = 'S1'><span style=' font-style: italic;'>A. Specifying Input</span><span> </span></h2><div  class = 'S2'><span>First, declare all the input variables such as file paths and filenames of input files and output files.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% File names and location</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >modelFileLocation = strcat(pwd, filesep); </span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >modelFileName = </span><span style="color: rgb(170, 4, 249);">'recon2_2_irr.mat'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >efmFileLocation = strcat(pwd, filesep);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >efmFileName = </span><span style="color: rgb(170, 4, 249);">'EFM-Dglucose-E10000-10757.txt'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >geneStatFlieLocation = strcat(pwd, filesep);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >geneStatFileName = </span><span style="color: rgb(170, 4, 249);">'data_useCase_2.txt'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Prefixes for the output files </span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >organism = </span><span style="color: rgb(170, 4, 249);">'human'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >dataFile = </span><span style="color: rgb(170, 4, 249);">'dataForMapping_EFM'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >efmBackboneFile = </span><span style="color: rgb(170, 4, 249);">'efmBackbone'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >efmPercentage = </span><span style="color: rgb(170, 4, 249);">'_EFMBackbone_percentage_'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >dataFileExt = </span><span style="color: rgb(170, 4, 249);">'.txt'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >sbmlExt = </span><span style="color: rgb(170, 4, 249);">'.xml'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >efmPre = </span><span style="color: rgb(170, 4, 249);">'_EFM_'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >metRem = </span><span style="color: rgb(170, 4, 249);">'_woUbMet_'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% Output file names</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >GSCFileName = strcat(</span><span style="color: rgb(170, 4, 249);">'GSC-'</span><span >, organism, dataFileExt); </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >GSSFileName = strcat(</span><span style="color: rgb(170, 4, 249);">'GSS-'</span><span >, organism, dataFileExt); </span></span></div></div></div><div  class = 'S9'><span>Also define the numeric indices of reactions of interest i.e., acetate release and glucose uptake, which will be used for calculating yield.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >roi_glcUp = 10757; </span><span style="color: rgb(2, 128, 9);">% glucose uptake (EX_glc(e)_b, also the active reaction for EFM Generation)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >roi_lacRel = 6613; </span><span style="color: rgb(2, 128, 9);">% lactate release, EX_lac_L(e)_f</span></span></div></div></div><div  class = 'S9'><span>If you know the reaction ID from the model, the corresponding index can be calculated as:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% index = find(contains(model_irr.rxns, 'EX_lac_L(e)_f')); % gives 6613 as output</span></span></div></div></div><div  class = 'S9'><span style=' font-weight: bold; font-style: italic;'>i. Import model, EFMs and gene level statistics</span></div><div  class = 'S2'><span>Load the same model that was used for generating set of EFMs. Here, Recon 2.2 was converted to a model containing only unidirectional (irreversible) reactions for EFM calculation.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >model_irr = importdata([modelFileLocation modelFileName]);</span></span></div></div></div><div  class = 'S9'><span>Alternatively, xml model file can also be used. In that case, use the COBRA function </span><span style=' font-style: italic;'>readCbModel.</span></div><div  class = 'S2'><span>Next, import the files containing EFMs. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >EFMRxns = efmImport(efmFileLocation, efmFileName);</span></span></div></div></div><div  class = 'S9'><span>Also, load the text file containing gene level statistics</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >geneData = readtable([geneStatFlieLocation geneStatFileName]); </span></span></div></div></div><div  class = 'S9'><span style=' font-weight: bold; text-decoration: underline;'>TIMING</span><span>: For large models and large set of EFMs, this step can take a few seconds to couple of minutes.</span></div><h2  class = 'S5'><span style=' font-style: italic;'>B. EFM Selection</span></h2><div  class = 'S2'><span style=' font-weight: bold; font-style: italic;'>i. Filter EFMs</span></div><div  class = 'S2'><span>Next, retain EFMs containing the acetate release reaction for further analysis.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >[filteredEFMs, filteredRows] = efmFilter(EFMRxns, roi_lacRel); </span></span></div></div></div><div  class = 'S9'><span>This results in 19 EFMs containing both glucose uptake and lactate release reactions. If large number of EFMs still remain, the above step can be repeated with additional reactions of interest.</span></div><div  class = 'S2'><span style=' font-weight: bold; font-style: italic;'>ii. EFM Enrichment</span></div><div  class = 'S2'><span>Enrichment is a two-step process consisting of preprocessing and enrichment analysis. Preprocessing will be done in MATLAB and enrichment itself will be done in R using the function 'runGSA()' which is a part of piano package [3]. Piano is useful for performing an enrichment analysis using custom gene sets.   </span></div><div  class = 'S2'><span style=' font-weight: bold; font-style: italic;'>a. Generate input files for EFM Enrichment (preprocessing)</span></div><div  class = 'S2'><span>The function runGSA() requires two input files - Gene Set Collection (GSC) and Gene Set Statistics (GSS). </span></div><div  class = 'S2'><span>GSC file: The gene set collection should describe the grouping of genes into gene sets. In this case, each EFM itself will be a set containing the reactions. A two-column (space-separated) file: &lt;EFM name or number&gt; and &lt;reaction ID&gt; will be generated.</span></div><div  class = 'S2'><span>GSS file:  The gene set collection should contain a p-value and expression value for each gene. Since our sets contain reactions, the gene level statistics will be mapped onto the reactions using the GPR rules defined in the model. A three-column (space-separated) file: &lt;reaction ID&gt; &lt;fold change&gt; &lt;p value&gt; will be generated.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% create the expreData variale with specific fields (gene, value, sig) as required by COBRA functions that map expression from genes to reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >exprData.gene = geneData.GeneID;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >exprData.value = geneData.log2FC;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >exprData.sig = geneData.PValue;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% exprSig = geneData.PValue;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >minSum = false;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >efmEnrichmentAnalysis1filteredEFMs, model_irr, exprData, GSCFileName, GSSFileName, minSum);</span></span></div></div></div><div  class = 'S9'><span>Once the input files needed for enrichment are written out, call the R script for performing enrichment.</span></div><div  class = 'S2'><span style=' font-weight: bold; font-style: italic;'>b. Call R script for Enrichment</span></div><div  class = 'S2'><span style=' font-weight: bold; text-decoration: underline;'>NOTE</span><span>: Modifiy the R directory according to the version of R installed on your system, if necessary.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% ! "C:\Program Files\R\R-3.5.1\bin\Rscript.exe" EFM_enrichment.R</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(178, 140, 0);">! "C:\Program Files\R\R-3.5.0\bin\Rscript.exe" EFM_enrichment.R</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="1FAF376D" data-scroll-top="null" data-scroll-left="null" data-width="428" data-height="227" data-hashorizontaloverflow="true" data-testid="output_0" style="max-height: 261px; width: 458px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Loading required package: piano 
Warning message:
package 'piano' was built under R version 3.5.1 
Running gene set analysis:
Checking arguments...done!
Final gene/gene-set association: 2192 genes and 19 gene sets
  Details:
  Calculating gene set statistics from 2192 out of 11567 gene-level statistics
  Using all 11567 gene-level statistics for significance estimation
  Removed 0 genes from GSC due to lack of matching gene statistics
  Removed 0 gene sets containing no genes after gene removal
  Removed additionally 0 gene sets not matching the size limits
  Loaded additional information for 0 gene sets
Calculating gene set statistics...done!
Calculating gene set significance...done!
Adjusting for multiple testing...done!</div></div></div></div></div><div  class = 'S9'><span>An output file "EFMEnrichmentResults.xls" will be written in the working directory. </span></div><div  class = 'S2'><span style=' font-weight: bold; text-decoration: underline;'>TIMING</span><span>: Enrichment - Preprocessing and running the R script could take a few minutes.</span></div><div  class = 'S2'><span style=' font-weight: bold; text-decoration: underline;'>TROUBLESHOOTING:</span><span> Check the location where R is installed in your system. If different from the code above, replace the path.</span></div><div  class = 'S2'><span style=' font-weight: bold; font-style: italic;'>iii. EFM Selection</span></div><div  class = 'S2'><span style=' font-weight: bold; font-style: italic;'>a. Get indices of enriched EFMs </span></div><div  class = 'S2'><span>Inspect the resulting file to select the most significant EFM(s).</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >upIndices = [83 29 42 28 51 84 41 69 70 79 82 62 36 35 12 39 64];</span></span></div></div></div><div  class = 'S9'><span style=' font-weight: bold; font-style: italic;'>b. Find Subsystem Occurrence in enriched EFMs and plot</span></div><div  class = 'S2'><span>Once you have identified the most significant EFMs, extract the distribution of subsystems to select an EFM to visualize.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >upEFMs = EFMRxns(upIndices, :);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >[subsysSummaryUp, uniqSubsysUp, countSubPerEFMUp] = efmSubsystemsExtraction(model_irr, upEFMs);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >tt=1:size(upEFMs, 1)</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >   len(tt) = length(find(upEFMs(tt,:))) ;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >countPerLenUp = countSubPerEFMUp ./ len;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >tt = 1:size(countPerLenUp,1)</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >    rowZeros(tt) = length(find(countPerLenUp(tt,:)==0));</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >countPerLen_filtered = countPerLenUp(rowZeros &lt; size(upEFMs,1)-round(0.5*size(upEFMs,1)),:);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span >uniqSubsysLeu_filtered = uniqSubsysUp(rowZeros &lt; size(upEFMs,1)-round(0.5*size(upEFMs,1)),:);</span></span></div></div></div><div  class = 'S2'><span style=' font-weight: bold; text-decoration: underline;'>CAUTION</span><span style=' text-decoration: underline;'>: DO NOT CLICK ON RUN! (FOR THIS SECTION ONLY)</span></div><div  class = 'S2'><span style=' text-decoration: underline;'>Select the below code for plotting and press F9</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >g3=clustergram(countPerLen_filtered);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >g3.Colormap=flipud(summer(500));</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >g3.Cluster=3;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >g3.Symmetric=false;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >g3.DisplayRange=60;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >g3.ColumnLabels=cellstr(strcat(</span><span style="color: rgb(170, 4, 249);">"EFM"</span><span >, string(upIndices)));</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >g3.RowLabels=cellstr(uniqSubsysLeu_filtered);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >g3.Colormap(1,:) = [1, 1, 1];</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >g3.ColumnLabelsRotate = 30;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >addTitle(g3, </span><span style="color: rgb(170, 4, 249);">"Subsystem Occurrence in EFMs containing Glucose uptake and Leucine release "</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div></div><div  class = 'S2'><span>As seen from the plot, EFM 83 is the shortest and stands out due to its highest density of reactions in Glycolysis/gluconeogenesis and Galactose metabolism and thus, was selected for visualization</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: pre"><span >selectedEFM = 83;</span></span></div></div></div><div  class = 'S9'><span>Once you have selected an EFM, you can extract the reactions in the EFM as a submodel (SMBL file) as explained below in C (Submodel Creation). This SBML file will be imported into Cytoscape in step D. You can choose to extract the SBML file with or without molecules that are connected to several reactions such as ATP and H2O using step C.i. or C.ii. respectively. We reccommend removing highly connected metabolites to improve the readability of the visualized network. In cases of analyzing/visualizing several EFMs, you can first visualize without ubquitous molecules to understand the underlying process and then look at the corresponding network with ubiquitous molecules.</span></div><div  class = 'S2'><span>An SBML file will be generated as output in each substep - C.i., C.ii. and C.iii. You can choose to run one step at a time. Follow the instructions in step D for each SBML file that you want to visualize.</span></div><h2  class = 'S1'><span style=' font-style: italic;'>C. Submodel Creation</span></h2><div  class = 'S2'><span style=' font-weight: bold; font-style: italic;'>i. Write submodel (with ubiquitous metabolites)</span></div><div  class = 'S2'><span>Get the indices of all reactions present in the selected EFM and extract it as a SBML file.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >selectedRxns = EFMRxns(selectedEFM, find(EFMRxns(selectedEFM,:))); </span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >efmSubmodelExtractionAsSBML(model_irr, selectedRxns, char(strcat(organism, efmPre, string(selectedEFM), sbmlExt)));</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="CFB54305" data-scroll-top="null" data-scroll-left="null" data-width="428" data-height="18" data-hashorizontaloverflow="false" data-testid="output_1" style="max-height: 261px; width: 458px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Document written</div></div></div></div></div><div  class = 'S9'><span>SBML file with all the metabolites is written in the current working directory.</span></div><div  class = 'S2'><span style=' font-weight: bold; text-decoration: underline;'>TIMING</span><span>: Writing SBML file could take a few seconds</span></div><div  class = 'S2'><span style=' font-weight: bold; font-style: italic;'>ii. Write submodel (withOUT ubiquitous metabolites)</span></div><div  class = 'S2'><span>Define the list of ubiquitous metabolites to be removed, making sure that the metabolite IDs match those in the model. Extract it as a SBML file</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >ubiquitousMets = {</span><span style="color: rgb(170, 4, 249);">'h'</span><span >;</span><span style="color: rgb(170, 4, 249);">'h2'</span><span >;</span><span style="color: rgb(170, 4, 249);">'h2o'</span><span >;</span><span style="color: rgb(170, 4, 249);">'atp'</span><span >;</span><span style="color: rgb(170, 4, 249);">'pi'</span><span >;</span><span style="color: rgb(170, 4, 249);">'adp'</span><span >;</span><span style="color: rgb(170, 4, 249);">'nad'</span><span >;</span><span style="color: rgb(170, 4, 249);">'ppi'</span><span >;</span><span style="color: rgb(170, 4, 249);">'pppi'</span><span >;</span><span style="color: rgb(170, 4, 249);">'nadh'</span><span >;</span><span style="color: rgb(170, 4, 249);">'nadp'</span><span >;</span><span style="color: rgb(170, 4, 249);">'nadph'</span><span >;</span><span style="color: rgb(170, 4, 249);">'amp'</span><span >;</span><span style="color: rgb(170, 4, 249);">'co2'</span><span >;</span><span style="color: rgb(170, 4, 249);">'coa'</span><span >;</span><span style="color: rgb(170, 4, 249);">'nh4'</span><span >;</span><span style="color: rgb(170, 4, 249);">'o2'</span><span >;</span><span style="color: rgb(170, 4, 249);">'fad'</span><span >;</span><span style="color: rgb(170, 4, 249);">'fadh2'</span><span >;</span><span style="color: rgb(170, 4, 249);">'fe2'</span><span >;</span><span style="color: rgb(170, 4, 249);">'fe3'</span><span >;</span><span style="color: rgb(170, 4, 249);">'na1'</span><span >;</span><span style="color: rgb(170, 4, 249);">'gdp'</span><span >;</span><span style="color: rgb(170, 4, 249);">'gtp'</span><span >;</span><span style="color: rgb(170, 4, 249);">'udp'</span><span >;</span><span style="color: rgb(170, 4, 249);">'ctp'</span><span >;</span><span style="color: rgb(170, 4, 249);">'utp'</span><span >;</span><span style="color: rgb(170, 4, 249);">'cdp'</span><span >;</span><span style="color: rgb(170, 4, 249);">'fmn'</span><span >;</span><span style="color: rgb(170, 4, 249);">'fmnh2'</span><span >};</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >efmSubmodelExtractionAsSBML(model_irr, selectedRxns, char(strcat(organism, efmPre, string(selectedEFM), metRem, sbmlExt)), 1, ubiquitousMets);</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="6451C618" data-scroll-top="null" data-scroll-left="null" data-width="428" data-height="18" data-hashorizontaloverflow="false" data-testid="output_2" style="max-height: 261px; width: 458px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Document written</div></div></div></div></div><div  class = 'S9'><span>SBML file withOUT all the specified metabolites is written in the current working directory.</span></div><div  class = 'S2'><span>In addition to looking at single EFMs, you can also use this workflow to extract an "EFM Backbone", i.e., commonly occurring reactions in the EFMs A certain percentage can be used as an input criteria , for example, 80%, 50%, 25%.</span></div><div  class = 'S2'><span style=' font-weight: bold; font-style: italic;'>iii. Extract reactions occurring in a certain percentage of EFMs</span></div><div  class = 'S2'><span>First, set the desired percentage and extract the reactions that are present in 'percentage' of EFMs. Input this set of reactions to the function 'efmSubmodelExtractionAsSBML' to write it out as an SBML file.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >percentage = 40;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >[selectedRxns,  rxnDistribution] = efmBackboneExtraction(filteredEFMs, percentage);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% x = h.BinEdges;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% y = h.Values;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% text(x(1:end-1)</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >rxnDistribution.rxnName = model_irr.rxns(rxnDistribution.rxnID);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S11'><span style="white-space: pre"><span >efmSubmodelExtractionAsSBML(model_irr, selectedRxns.rxnID, char(strcat(organism, efmPercentage, string(percentage), sbmlExt)), 1, ubiquitousMets);</span></span></div><div  class = 'S12'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="97271A18" data-scroll-top="null" data-scroll-left="null" data-width="428" data-height="18" data-hashorizontaloverflow="false" data-testid="output_3" style="max-height: 261px; width: 458px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Document written</div></div></div></div><div class="inlineWrapper"><div  class = 'S13'></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% % [OPTIONAL] Plot the number of reactions present in different proportion of EFMs and save this data in a text and map on visualized network</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% figure; h = his,y,num2str(y'),'vert','bottom','horiz','left');</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% writetable(rxnDistribution, strcat(organism, dataFile, '_', efmBackboneFile, dataFileExt));</span></span></div></div></div><div  class = 'S9'><span>Note: The tutorial on use case 1 included an additional optional step here to extract flux data for mapping. Apart from flux data, you can also map gene expression data on the visualized network for which, the same input file from step B can be used.</span><span style=' font-weight: bold; font-style: italic;'> </span></div><h2  class = 'S5'><span>D. Visualization</span></h2><div  class = 'S2'><span>Follow these steps for each submodel SBML file that you wish to visualize.</span></div><ul  class = 'S3'><li  class = 'S4'><span>Make sure Cytoscape (v3 and above) is installed along with cy3sbml (app) and yFiles layout</span></li><li  class = 'S4'><span>Open the SBML file using the cy3sbml app. [Click on the cy3sbl app-specific icon and select the EFM submodel file to open. For every file imported, three networks are generated and are identified using the prefixes: All, Base and Kinetic. Use the Base network for rest of the tutorial.</span></li></ul><div  class = 'S2'><img class = "imageNode" src = "" alt = "" style = "vertical-align: baseline"></img></div><ul  class = 'S3'><li  class = 'S4'><span>Import the style file supplied with this workflow (efmviz_style.xml) and apply this style on the Base network. You can now see a network with three different types of nodes representing reactions, genes and metabolites and edges representing interaction between these molecule types. This is how your EFM looks!</span></li><li  class = 'S4'><span>Set the layout using Layout -&gt; yFiles Orthogonal layout. In our experience orthogonal layout works best for EFMs with 50-60 reactions. Beyond that, the EFM can be visualized, but will be difficult to interpret.</span></li><li  class = 'S4'><span>(Optionally) Import the data file generated in step C.iv to map on the reaction edges. Use the R script (see below for details)</span></li></ul><div  class = 'S2'><span>For running this step on several files, follow instructions in the R script, efmviz_rcy3_human.Rmd. </span></div><div  class = 'S2'><span style=' font-weight: bold; text-decoration: underline;'>Try it yourself</span><span style=' font-weight: bold;'>:</span><span> Repeat steps C and D for other EFMs.</span></div><h2  class = 'S1'><span>REFERENCES</span></h2><div  class = 'S2'><span>[1] Feist, A.M.; Henry, C.S.; Reed, J.L.; Krummenacker, M.; Joyce, A.R.; Karp, P.D.; Broadbelt, L.J.;Hatzimanikatis, V.; Palsson, B.Ø. A genome-scale metabolic reconstruction for Escherichia coli K-12 MG1655 that accounts for 1260 ORFs and thermodynamic information. Molecular Systems Biology 2007, 3, 121. doi:10.1038/MSB4100155.</span></div><div  class = 'S2'><span>[2] Swainston, N.; Smallbone, K.; Hefzi, H.; Dobson, P.D.; Brewer, J.; Hanscho, M.; Zielinski, D.C.; Ang, K.S.; Gardiner, N.J.; Gutierrez, J.M.; Kyriakopoulos, S.; Lakshmanan, M.; Li, S.; Liu, J.K.; Martínez, V.S.; Orellana, C.A.; Quek, L.E.; Thomas, A.; Zanghellini, J.; Borth, N.; Lee, D.Y.; Nielsen, L.K.; Kell, D.B.; Lewis, N.E.; Mendes, P. Recon 2.2: from reconstruction to model of human metabolism. Metabolomics 2016, 12, 109. doi:10.1007/s11306-016-1051-4.</span></div><div  class = 'S2'><span>[3] Varemo, L., Nielsen, J., Nookaew, I. Enriching the gene set analysis of genomewide data by incorporating directionality of gene expression and combining statistical hypotheses and methods. Nucleic Acids Research 2013, 41, 8. doi:10.1093/nar/gkt111</span></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% Visualize Elementary Flux Modes from Human model 
%% Author(s): Chaitra Sarathy (Maastricht University, The Netherlands)
% 
%% Reviewer(s): 
% _Include those who have reviewed the tutorial, and are therefore partly responsible 
% for picking up errors or omissions._
%% INTRODUCTION
% Elementary Flux Modes (EFMs) are a tool for constraint-based modelling and 
% metabolic network analysis. However, systematic and automated visualization 
% of EFMs, capable of integrating various data types is still a challenge. This 
% tutorial is part of the study by Sarathy et al. (2019, in preparation), "EFMviz: 
% A COBRA Toolbox extension to visualize Elementary Flux Modes in genome-scale 
% metabolic models". 
% 
% This workflow can be used to visualize EFMs calculated using any currently 
% available methods. Inputs for the workflow are: (1) the text file containing 
% EFMs and (2) the genome-scale metabolic model used for generating EFMs. Additionally, 
% omics data can be used to select one (or more) EFM(s) for visualization among 
% the vast number of EFMs that usually results from EFM calculation. Sarathy et 
% al. (2019, in preparation) demonstrate two examples of such a selection procedure 
% (a) Yield Analysis and (b) EFM Enrichment with EFMs calculated from genome-scale 
% metabolic models, (1) iAF1260 (E. coli) [1] and (2) Recon 2.2 (Human) [2]. 
% 
% Each use case is explained with a tutorial. This tutorial corresponds to the 
% second use case and it describes how to analyze and visualize EFMs from the 
% human model Recon 2.2_._ The aim is to apply the workflow on a large-scale model 
% and it can be used as a means to study cancer metabolism. EFMs were generated 
% for glucose uptake reaction. For details on how EFMs were generated, please 
% refer the citation section at <https://github.com/macsbio/efmviz https://github.com/macsbio/efmviz>
% 
% The steps in the workflow include:
% 
% *A. Specifying Inputs*
% 
% i. Import input data: file containing EFMs (.txt), model file (.mat), file 
% containing gene level statistics typically resulting from a differential gene 
% expression analysis (.txt)
% 
% *B. EFM Selection*
% 
% i. Filter EFMs containing a reaction of interest
% 
% ii. Perform EFM Enrichment Analysis
% 
% a. Generate input files for EFM Enrichment (preprocessing)
% 
% b. Call R script for Enrichment (EFM_enrichment.R)
% 
% iii. EFM Selection
% 
% a. Get indices of enriched EFMs 
% 
% b. Find Subsystem Occurrence in enriched EFMs for selecting an EFM for visualization, 
% for example, EFM with reactions in TCA Cycle
% 
% *C. Submodel Creation*
% 
% i. Write submodel (with ubiquitous metabolites)
% 
% ii. Write submodel (withOUT ubiquitous metabolites)
% 
% iii. Extract reactions occurring in a certain percentage of EFMs
% 
% *D. Visualization*
% 
% Open the SBML file in Cytoscape and map the data used during step B.
% 
% The following sections will walk you through the code for performing the steps 
% A-D, thereby reproducing the main results from Sarathy et al (2019, in preparation) 
% (Use Case 2). For more details on the workflow please refer <https://github.com/macsbio/efmviz 
% https://github.com/macsbio/efmviz>
% 
% *Input File Formats*
%% 
% * model (.mat): same as the model used for EFM generation 
% * file containing EFMs (.txt): each row corresponds to an EFM and column has 
% indices of reactions active in the EFM. Each row can have varied number of columns. 
% The first column should contain an EFM ID.
% * file containing gene level statistics (.txt): Tab separated file containing 
% an expression level (log2 fold change) and p-value for every gene. Columns: 
% <GeneID><GeneName><log2FC><PValue>. Make sure to have the gene IDs in the same 
% format as in the model.
%% PROCEDURE
%% _A. Specifying Input_ 
% First, declare all the input variables such as file paths and filenames of 
% input files and output files.

% File names and location
modelFileLocation = strcat(pwd, filesep); 
modelFileName = 'recon2_2_irr.mat';

efmFileLocation = strcat(pwd, filesep);
efmFileName = 'EFM-Dglucose-E10000-10757.txt';

geneStatFlieLocation = strcat(pwd, filesep);
geneStatFileName = 'data_useCase_2.txt';

% Prefixes for the output files 
organism = 'human';
dataFile = 'dataForMapping_EFM';
efmBackboneFile = 'efmBackbone';
efmPercentage = '_EFMBackbone_percentage_';
dataFileExt = '.txt';
sbmlExt = '.xml';
efmPre = '_EFM_';
metRem = '_woUbMet_';

% Output file names
GSCFileName = strcat('GSC-', organism, dataFileExt); 
GSSFileName = strcat('GSS-', organism, dataFileExt); 
%% 
% Also define the numeric indices of reactions of interest i.e., acetate release 
% and glucose uptake, which will be used for calculating yield.

roi_glcUp = 10757; % glucose uptake (EX_glc(e)_b, also the active reaction for EFM Generation)
roi_lacRel = 6613; % lactate release, EX_lac_L(e)_f
%% 
% If you know the reaction ID from the model, the corresponding index can be 
% calculated as:

% index = find(contains(model_irr.rxns, 'EX_lac_L(e)_f')); % gives 6613 as output
%% 
% _*i. Import model, EFMs and gene level statistics*_
% 
% Load the same model that was used for generating set of EFMs. Here, Recon 
% 2.2 was converted to a model containing only unidirectional (irreversible) reactions 
% for EFM calculation.

model_irr = importdata([modelFileLocation modelFileName]);
%% 
% Alternatively, xml model file can also be used. In that case, use the COBRA 
% function _readCbModel._
% 
% Next, import the files containing EFMs. 

EFMRxns = efmImport(efmFileLocation, efmFileName);
%% 
% Also, load the text file containing gene level statistics

geneData = readtable([geneStatFlieLocation geneStatFileName]); 
%% 
% *TIMING*: For large models and large set of EFMs, this step can take a few 
% seconds to couple of minutes.
%% _B. EFM Selection_
% _*i. Filter EFMs*_
% 
% Next, retain EFMs containing the acetate release reaction for further analysis.

[filteredEFMs, filteredRows] = efmFilter(EFMRxns, roi_lacRel); 
%% 
% This results in 19 EFMs containing both glucose uptake and lactate release 
% reactions. If large number of EFMs still remain, the above step can be repeated 
% with additional reactions of interest.
% 
% _*ii. EFM Enrichment*_
% 
% Enrichment is a two-step process consisting of preprocessing and enrichment 
% analysis. Preprocessing will be done in MATLAB and enrichment itself will be 
% done in R using the function 'runGSA()' which is a part of piano package [3]. 
% Piano is useful for performing an enrichment analysis using custom gene sets.   
% 
% _*a. Generate input files for EFM Enrichment (preprocessing)*_
% 
% The function runGSA() requires two input files - Gene Set Collection (GSC) 
% and Gene Set Statistics (GSS). 
% 
% GSC file: The gene set collection should describe the grouping of genes into 
% gene sets. In this case, each EFM itself will be a set containing the reactions. 
% A two-column (space-separated) file: <EFM name or number> and <reaction ID> 
% will be generated.
% 
% GSS file:  The gene set collection should contain a p-value and expression 
% value for each gene. Since our sets contain reactions, the gene level statistics 
% will be mapped onto the reactions using the GPR rules defined in the model. 
% A three-column (space-separated) file: <reaction ID> <fold change> <p value> 
% will be generated.

% create the expreData variale with specific fields (gene, value, sig) as required by COBRA functions that map expression from genes to reactions
exprData.gene = geneData.GeneID;
exprData.value = geneData.log2FC;
exprData.sig = geneData.PValue;
% exprSig = geneData.PValue;
minSum = false;

efmEnrichmentAnalysis1filteredEFMs, model_irr, exprData, GSCFileName, GSSFileName, minSum);
%% 
% Once the input files needed for enrichment are written out, call the R script 
% for performing enrichment.
% 
% _*b. Call R script for Enrichment*_
% 
% *NOTE*: Modifiy the R directory according to the version of R installed on 
% your system, if necessary.

% ! "C:\Program Files\R\R-3.5.1\bin\Rscript.exe" EFM_enrichment.R
! "C:\Program Files\R\R-3.5.0\bin\Rscript.exe" EFM_enrichment.R
%% 
% An output file "EFMEnrichmentResults.xls" will be written in the working directory. 
% 
% *TIMING*: Enrichment - Preprocessing and running the R script could take a 
% few minutes.
% 
% *TROUBLESHOOTING:* Check the location where R is installed in your system. 
% If different from the code above, replace the path.
%% 
% _*iii. EFM Selection*_
% 
% _*a. Get indices of enriched EFMs*_ 
% 
% Inspect the resulting file to select the most significant EFM(s).

upIndices = [83 29 42 28 51 84 41 69 70 79 82 62 36 35 12 39 64];
%% 
% _*b. Find Subsystem Occurrence in enriched EFMs and plot*_
% 
% Once you have identified the most significant EFMs, extract the distribution 
% of subsystems to select an EFM to visualize.

upEFMs = EFMRxns(upIndices, :);
[subsysSummaryUp, uniqSubsysUp, countSubPerEFMUp] = efmSubsystemsExtraction(model_irr, upEFMs);

for tt=1:size(upEFMs, 1)
   len(tt) = length(find(upEFMs(tt,:))) ;
end

countPerLenUp = countSubPerEFMUp ./ len;
for tt = 1:size(countPerLenUp,1)
    rowZeros(tt) = length(find(countPerLenUp(tt,:)==0));
end
countPerLen_filtered = countPerLenUp(rowZeros < size(upEFMs,1)-round(0.5*size(upEFMs,1)),:);
uniqSubsysLeu_filtered = uniqSubsysUp(rowZeros < size(upEFMs,1)-round(0.5*size(upEFMs,1)),:);
%% 
% *CAUTION*: DO NOT CLICK ON RUN! (FOR THIS SECTION ONLY)
% 
% Select the below code for plotting and press F9

g3=clustergram(countPerLen_filtered);
g3.Colormap=flipud(summer(500));
g3.Cluster=3;
g3.Symmetric=false;
g3.DisplayRange=60;
g3.ColumnLabels=cellstr(strcat("EFM", string(upIndices)));
g3.RowLabels=cellstr(uniqSubsysLeu_filtered);
g3.Colormap(1,:) = [1, 1, 1];
g3.ColumnLabelsRotate = 30;
addTitle(g3, "Subsystem Occurrence in EFMs containing Glucose uptake and Leucine release ");

%% 
% As seen from the plot, EFM 83 is the shortest and stands out due to its highest 
% density of reactions in Glycolysis/gluconeogenesis and Galactose metabolism 
% and thus, was selected for visualization

selectedEFM = 83;
%% 
% Once you have selected an EFM, you can extract the reactions in the EFM as 
% a submodel (SMBL file) as explained below in C (Submodel Creation). This SBML 
% file will be imported into Cytoscape in step D. You can choose to extract the 
% SBML file with or without molecules that are connected to several reactions 
% such as ATP and H2O using step C.i. or C.ii. respectively. We reccommend removing 
% highly connected metabolites to improve the readability of the visualized network. 
% In cases of analyzing/visualizing several EFMs, you can first visualize without 
% ubquitous molecules to understand the underlying process and then look at the 
% corresponding network with ubiquitous molecules.
% 
% An SBML file will be generated as output in each substep - C.i., C.ii. and 
% C.iii. You can choose to run one step at a time. Follow the instructions in 
% step D for each SBML file that you want to visualize.
%% _C. Submodel Creation_
% _*i. Write submodel (with ubiquitous metabolites)*_
% 
% Get the indices of all reactions present in the selected EFM and extract it 
% as a SBML file.

selectedRxns = EFMRxns(selectedEFM, find(EFMRxns(selectedEFM,:))); 
efmSubmodelExtractionAsSBML(model_irr, selectedRxns, char(strcat(organism, efmPre, string(selectedEFM), sbmlExt)));
%% 
% SBML file with all the metabolites is written in the current working directory.
% 
% *TIMING*: Writing SBML file could take a few seconds
% 
% _*ii. Write submodel (withOUT ubiquitous metabolites)*_
% 
% Define the list of ubiquitous metabolites to be removed, making sure that 
% the metabolite IDs match those in the model. Extract it as a SBML file

ubiquitousMets = {'h';'h2';'h2o';'atp';'pi';'adp';'nad';'ppi';'pppi';'nadh';'nadp';'nadph';'amp';'co2';'coa';'nh4';'o2';'fad';'fadh2';'fe2';'fe3';'na1';'gdp';'gtp';'udp';'ctp';'utp';'cdp';'fmn';'fmnh2'};
efmSubmodelExtractionAsSBML(model_irr, selectedRxns, char(strcat(organism, efmPre, string(selectedEFM), metRem, sbmlExt)), 1, ubiquitousMets);
%% 
% SBML file withOUT all the specified metabolites is written in the current 
% working directory.
% 
% In addition to looking at single EFMs, you can also use this workflow to extract 
% an "EFM Backbone", i.e., commonly occurring reactions in the EFMs A certain 
% percentage can be used as an input criteria , for example, 80%, 50%, 25%.
% 
% _*iii. Extract reactions occurring in a certain percentage of EFMs*_
% 
% First, set the desired percentage and extract the reactions that are present 
% in 'percentage' of EFMs. Input this set of reactions to the function 'efmSubmodelExtractionAsSBML' 
% to write it out as an SBML file.

percentage = 40;
[selectedRxns,  rxnDistribution] = efmBackboneExtraction(filteredEFMs, percentage);
% x = h.BinEdges;
% y = h.Values;
% text(x(1:end-1)
rxnDistribution.rxnName = model_irr.rxns(rxnDistribution.rxnID);
efmSubmodelExtractionAsSBML(model_irr, selectedRxns.rxnID, char(strcat(organism, efmPercentage, string(percentage), sbmlExt)), 1, ubiquitousMets);

% % [OPTIONAL] Plot the number of reactions present in different proportion of EFMs and save this data in a text and map on visualized network
% figure; h = his,y,num2str(y'),'vert','bottom','horiz','left');
% writetable(rxnDistribution, strcat(organism, dataFile, '_', efmBackboneFile, dataFileExt));
%% 
% Note: The tutorial on use case 1 included an additional optional step here 
% to extract flux data for mapping. Apart from flux data, you can also map gene 
% expression data on the visualized network for which, the same input file from 
% step B can be used. 
%% D. Visualization
% Follow these steps for each submodel SBML file that you wish to visualize.
%% 
% * Make sure Cytoscape (v3 and above) is installed along with cy3sbml (app) 
% and yFiles layout
% * Open the SBML file using the cy3sbml app. [Click on the cy3sbl app-specific 
% icon and select the EFM submodel file to open. For every file imported, three 
% networks are generated and are identified using the prefixes: All, Base and 
% Kinetic. Use the Base network for rest of the tutorial.
%% 
% 
%% 
% * Import the style file supplied with this workflow (efmviz_style.xml) and 
% apply this style on the Base network. You can now see a network with three different 
% types of nodes representing reactions, genes and metabolites and edges representing 
% interaction between these molecule types. This is how your EFM looks!
% * Set the layout using Layout -> yFiles Orthogonal layout. In our experience 
% orthogonal layout works best for EFMs with 50-60 reactions. Beyond that, the 
% EFM can be visualized, but will be difficult to interpret.
% * (Optionally) Import the data file generated in step C.iv to map on the reaction 
% edges. Use the R script (see below for details)
%% 
% For running this step on several files, follow instructions in the R script, 
% efmviz_rcy3_human.Rmd. 
% 
% *Try it yourself:* Repeat steps C and D for other EFMs.
%% REFERENCES
% [1] Feist, A.M.; Henry, C.S.; Reed, J.L.; Krummenacker, M.; Joyce, A.R.; Karp, 
% P.D.; Broadbelt, L.J.;Hatzimanikatis, V.; Palsson, B.Ø. A genome-scale metabolic 
% reconstruction for Escherichia coli K-12 MG1655 that accounts for 1260 ORFs 
% and thermodynamic information. Molecular Systems Biology 2007, 3, 121. doi:10.1038/MSB4100155.
% 
% [2] Swainston, N.; Smallbone, K.; Hefzi, H.; Dobson, P.D.; Brewer, J.; Hanscho, 
% M.; Zielinski, D.C.; Ang, K.S.; Gardiner, N.J.; Gutierrez, J.M.; Kyriakopoulos, 
% S.; Lakshmanan, M.; Li, S.; Liu, J.K.; Martínez, V.S.; Orellana, C.A.; Quek, 
% L.E.; Thomas, A.; Zanghellini, J.; Borth, N.; Lee, D.Y.; Nielsen, L.K.; Kell, 
% D.B.; Lewis, N.E.; Mendes, P. Recon 2.2: from reconstruction to model of human 
% metabolism. Metabolomics 2016, 12, 109. doi:10.1007/s11306-016-1051-4.
% 
% [3] Varemo, L., Nielsen, J., Nookaew, I. Enriching the gene set analysis of 
% genomewide data by incorporating directionality of gene expression and combining 
% statistical hypotheses and methods. Nucleic Acids Research 2013, 41, 8. doi:10.1093/nar/gkt111
##### SOURCE END #####
-->
</div></body></html>